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Abstract 

In this paper, for the first time a theory is formulated that predicts 
velocity and spatial correlations between occupation numbers that oc- 
cur in lattice gas automata violating semi-detailed balance. Starting 
from a coupled BBGKY hierarchy for the n-particle distribution func- 
tions, cluster expansion techniques are used to derive approximate 
kinetic equations. In zeroth approximation the standard nonlinear 
Boltzmann equation is obtained; the next approximation yields the 
ring kinetic equation, similar to that for hard sphere systems, describ- 
ing the time evolution of pair correlations. As a quantitative test we 
calculate equal time correlation functions in equilibrium for two mod- 
els that violate semi-detailed balance. One is a model of interacting 
random walkers on a line, the other one is a two-dimensional fluid 
type model on a triangular lattice. The numerical predictions agree 
very well with computer simulations. 

Key words: non-Gibbs states, lack of detailed balance, static pair 
correlations, lattice gas automata, BBGKY hierarchy 
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1 Introduction 



In the theory of lattice gas automata (LGA's) the Gibbs distribution is known 
to be the unique equilibrium distribution if the dynamics satisfy the so-called 
semi-detailed balance ]I[ or Stueckelberg condition. The Gibbs distribu- 
tion depends on the phase space variables only through globally conserved 
quantities, such as the total number of particle N, the total momentum P 
or for thermal models the total energy H. If the collision rules are strictly 
local, involving particles on the same node only, then the Gibbs equilibrium 
state is completely factorized, i.e. the occupations of velocity channels at the 
same or at different nodes are completely uncorrelated. 

However, in LGA's that violate the condition of semi-detailed balance 
equal time correlations are know to exists in equilibrium between occupa- 
tion numbers of different velocity channels on the same or on different nodes 
[@- 0! Hi • The existence of such equilibrium correlations prevents the equi- 
librium distribution from being the Gibbs distribution, and therefore stan- 
dard equilibrium statistical mechanics does not apply to LGA's that violate 
(semi)-detailed balance. In the existing literature a quantitative or qualita- 
tive understanding of the structure of the equilibrium distribution for such 
models is entirely lacking. Fundamental questions are: is there a unique equi- 
librium state, or can there be several stationary states? To what extent does 
the final state depend on the initial state and on the details of the collision 
rules? How does the system approach the correlated equilibrium state. How 
can one calculate transport properties and time correlation functions? 

The aim of the present paper is to provide a theory for calculating time 
dependent distribution and correlation functions in LGA's that violate semi- 
detailed balance. Using the general ideas developed in the classical theories of 
fluids, we construct a BBGKY hierarchy for the distribution functions, and 
obtain by means of cluster expansion techniques a generalized Boltzmann 
equation for the single particle distribution function, coupled to the so-called 
ring equation for the pair correlation function. We are specially interested in 
calculating the structure of a spatially uniform equilibrium state. 

The organization of the paper is as follows. In section 2 we derive the 
first two equations of a coupled hierarchy for the joint distribution functions, 
and we introduce a cluster expansion around mean field theory that leads to 
a closure of the hierarchy at the level of two-point correlation functions. In 
section 3 the coupled kinetic equations are specialized to spatially uniform 
stationary states, and equations are constructed for the stationary values of 
the single particle distribution function and the pair correlation function; an 
iterative scheme for evaluating these equations numerically is discussed in 
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section 4. To establish the accuracy of the theory, we apply in section 5 the 
theory to a model of interacting random walkers on a line, and to a fluid- 
type model on the triangular lattice, and compare the theoretical predictions 
with detailed computer simulations performed for both models. We end with 
some conclusions and remarks about future developments in section 6. 

2 Kinetic theory 

2.1 Microdynamic equation 

We consider a lattice gas automaton defined on a regular d- dimensional lat- 
tice £ with periodic boundary conditions, containing V = L d nodes. On 
each node r e C there exist b allowed velocity channels Ci, c 2 , . . . , c&, cor- 
responding to the nearest neighbor lattice vectors; there may also be rest 

particles with |q| = 0. The system evolves at discrete time steps t — 0, 1, 

The microscopic configuration of the system at time t is given in terms of 
occupation numbers rii(r,t) = {0,1}, denoting the absence or presence of 
a particle in velocity channel (r, Cj). The state of node r is denoted by 
n(r,t) = {m(r,t);i = 1,2, ..,&}. 

One time step in the evolution involves an instantaneous collision step 
which transforms a precollision state n(r, t) into a postcollision state n*(r, t) 
on all nodes r independently. It is followed by a propagation step, 

ni (r + Cl ,t+l)=n*(r,t), (2.1) 

during which a particle at node r is moved to the nearest neighbor node 
r + Ci in the direction of its velocity Cj. The collision step is specified in 
terms of a matrix of transition probabilities A sa from a precollision or in- 
state s = {si, S2, • • • , Sb} to a postcollision or out-state a = {o"i, a 2 , ■ ■ ■ , o~b}- 
The matrix A sa satisfies the normalization condition 

£A S(T = 1. (2.2) 

cr 

There are additional constraints imposed on the transition matrix A sa due 
to the presence of local conservation laws, which are most easily expressed 
in terms of the collisional invariants a« = {1, q, |c?, ...}, as 

£ OiUiir, t) = a%n*(r, t). (2.3) 

i i 

The local conservation laws imply global conservation of the total quantities 
H = {N, P,H, . . .} = J2r,i CLiTii(r, t). In diffusive LGA's only the total num- 
ber of particles, N, is conserved; in fluid-type LGA's the total momentum, 



3 



P, is conserved in addition. There also exist thermal models where the total 
energy H is conserved. Furthermore, the majority of LGA's have additional 
(spurious) global invariants which are mostly 'staggered' in space and 
time. They can be included in the present discussion in a similar manner. 

Much is known about the properties of LGA's that satisfy the condition 
of semi- detailed balance 

£A SCT = 1, (2.4) 

s 

also known as the Stueckelberg condition [[J, or the stronger detailed balance 
condition, 

A S(7 = A as . (2-5) 

which implies the semi-detailed balance condition ( |2.4| ) on account of the 
normalization condition (|2.2|) . If condition (|2.4j) is satisfied then the equilib- 
rium distribution for the model is the Gibbs distribution, which only depends 
on the globally conserved quantities 7i.[| 

A comment on the difference between semi-detailed balance and detailed 
balance is appropriate. We group all states s(r) that can be mapped into one 
another by lattice symmetries (rotations, reflections) into the same equiva- 
lence class [§J. Then, suppose that for any given in-state, s(r), transitions 
are allowed to out-states, cr(r), belonging to at most M different equivalence 
classes. As long as M < 2 there is no distinction between semi-detailed 
balance and detailed balance. Almost all LGA's used in the literature have 
M < 2. Notable exceptions are the FCHC models §. 

In what follows, no condition apart from normalization will be imposed 
on the transition matrix A scr . Of course, all the results that will be obtained 
in the rest of this paper will hold for LGA's that satisfy semi-detailed balance 
as well, being just a special case of the more general class of LGA's violating 
semi-detailed balance. 



2.2 Hierarchy equations 

To construct the equations of motion for the occupation numbers we observe 
that the occupation numbers n(r,t) and n*(r,t) are related through the 
transition matrix A sa as, 

K(r,t) = J2<nA a J(sMr,t)), (2.6) 

as 

1 The semi-detailed balance and detailed balance conditions with respect to an arbitrary 
phase space distribution function Pq(s) have respectively the forms, J2 S Po($)A sa = Po{a) 
and Po(s)A sa = A as P^{a). The phase space distribution in the Gibbs state depends on 
the variables s(r) only through locally conserved quantities. Consequently Pq{s) — Pq{u) 
and the (semi)-detailed balance conditions reduce to the equations in the text. 
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where 8(s,n) = Ylj 8( s j, n j) is a product of b Kronecker delta functions. For 
the product of two occupation numbers on the same node we have 

n*{r,t)n*(r,t) =^2<r i <T j A aa S{s,n{r,t)). (2.7) 

as 

For every node r and every time t the Boolean realization A sa = A sa (r, t) = 
{0, 1} of the (s — > a)-transition is drawn from an (r, /^-independent proba- 
bility distribution with expectation value ((A sa )) = A sa . 

The hierarchy equations describe the time dependence of the distribution 
functions, i.e. the expectation value of products of occupation numbers. The 
single particle and two particle distribution functions, 



fi(r,t) = {m(r,t)) 
m r'. I) = (m(r,t)n v 



ft?\r,r',t) = (n i (r,t)n J -(r , ,f)), (2.8) 



are defined as averages (...) over an arbitrary distribution of initial occupa- 
tion numbers, {n(r, 0)}. Note that (r,r,t) = fi(r,t). 

We now introduce the single node and two node distribution functions, 

p(s,r,t) = (5(s,n(r,t))) (2.9) 
p^(s,r,s',r',t) = (S(s,n(r,t))5(s',n(r',t))), (2.10) 

By taking the average of Q2.6|) over an arbitrary distribution of initial states 
{n,(r, 0)}, as well as over the A-distribution of realizations of the transition 
matrix, and using ([2.1Q, we obtain 

fi(r + a,t + l) = (n*(r,t)) = J2^iAsMs,r,t). (2.11) 

as 

To obtain an analogous equation for the pair distribution function we must 
distinguish between the cases r — r' and r ^ r'. When r ^ r' the collisions 
for channels (r, Cj) and (r',Cj) are independent, and given by (|2.6| ). When 
r = r' however we must use (12.71). The combined result can be written as 



f%\r + Cl ,r + c v t + l) = (n*(r,t)n*(r',t)) 



[l-8(r,r')]l Yl ^A sa A sl(7/ p^(s,r,s',r',t)\ (2.12) 

^ as,a's' 

+ 5(r,r') |53<r i (7 i ,A«TP( s : r »*)| • 
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The above equations Q2.11Q and Q2.12|) constitute respectively the first two 
equations of an open hierarchy of coupled equations for the distribution func- 
tions, similar to the BBGKY hierarchy for continuous systems. It is straight- 
forward to extend this procedure to construct the hierarchy equations for 
with £ > 3, although the complexity increases rapidly. In this paper the first 
two equations of the BBGKY hierarchy suffice for our purpose. 

The node distribution functions are related to the ^-particle distribution 
functions in a simple manner. This follows from the following identity, valid 
for binary variables, 

S(s,n) = n*(«i,"i) = IK'C 1 " %) 1-SJ - ( 2 - 13 ) 

j 3 

As Sj takes only values or 1, the right hand side is a sum of products of at 
most b occupation numbers. The average of ( |2.13| ), p(s,r,t), is therefore a 
linear combination of ^-particle distribution functions /™ with £ = 1, 2, .., b, 
all referring to the same node, and with coefficients that are either +1 or — 1. 
Similarly is a linear combination of distribution functions. 

To understand the differences between the BBGKY hierarchy for dis- 
crete LGA's and continuous systems we note the following. In continuous 
systems usually only two-body interactions (additive forces) are considered; 
consequently the RHS of the £-th hierarchy equation contains only (£ + 1)- 
particle distribution functions. If one would include s-body interactions 
(s = 2, 3, 4, .., b) the £-th hierarchy equation would contain (£ + b— l)-particle 
distribution functions. In the LGA case 6-body interactions simultaneously 
occur on each node r, where b is the number of velocity channels at a node. 
Therefore the first hierarchy equation already involves up to 6-particle dis- 
tribution functions; the second equation even involves up to 26-particle func- 
tions. 

2.3 Cluster expansion 

An approximate closure of the hierarchy equations can be obtained by making 
a cluster expansion in terms of two-, three-, etc. point correlation functions 
and retaining correlation functions up to a certain order while neglecting the 
higher order ones. This procedure leads to approximate kinetic equations. 

The cluster functions G are defined through the so-called Ursell expan- 
sion, 

f§\r,r>) = mfjW+G^r') 

f$(r,r',r") = / 4 (r)/.,(rU(0 + fi(r)G j£ (r' ,r") + f 3 (r')G ie (r,r") 

+ fi{r")G i:i (r, r') + G ijt (r, r', r"), (2.14) 
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etc, where in general all channels (r, q), (r', Cj), and (r", Cj) are different. 
Solution of the recursion relations yields the pair and triplet functions in 
terms of the f^'. The cluster functions conveniently can be written in terms 
of fluctuations, 

5ni(r,t) = n l (r,t) - fi{r,t). (2.15) 
For instance, for the two point function we have 

G tJ (r,r',t) = fV(r,r',t) - / i (r,t)/ i (r',t) = (MM<^V, *)>, (2.16) 

and similarly we can write GV^(r,r' ,r" ,t) = (5ni(r,t)5nj(r',t)5ne(r",t)). 
It is important to note that the diagonal elements for (r, c$) = (r', Cj) are 
completely determined by fi(r,t), due to the Boolean character of rii(r,t), 

G«(r,r,f) = ((Mr,t)) 2 > = /i(r,*)(l - £(»",*)) = <*(»",*). (2.17) 

To perform the cluster expansion of the hierarchy equations the cluster 
expansion of the node distribution functions p(s) and p^ 2 \s,s') is needed. 
In appendix A the expansion up to terms linear in the pair correlations is 
worked out detail. It amounts to obtaining approximate expressions for p 
and p^ in terms of / and G. For the single node distribution function the 
results is 

p(s, r, t) = F(s, r, t) (l + £ ^f^ G^r, r,t) + ..\, (2.18) 

{ $<i 9k{r,t)g l {r,t) J 

as shown in appendix A. Here F(s,r,t) is the completely factorized single 
node distribution function, 

F(s,r,t) = n(/i(^)) Si (l - .Dir.')) 1 (2.19) 



In a similar way we can expand the two-node distribution function (see 
eq. (|A.8|) of appendix A). Substitution of these expansions in the hierarchy 
equations yields coupled approximate equations of motion for fi(r,t) and 
GV,(r,r',t), given by (|A~5| ) together with Q , (|A~T0|) , (|A~lp and (M - 
The kinetic equation in zeroth approximation is obtained by neglecting 
Gjj in ( |2.18| ), or equivalently by replacing the distribution functions p(s, r, t) 
in ( p. 11] ) by the fully factorized F(s,r,t). In that case the first hierarchy 
equation reduces to the nonlinear Boltzmann equation, 

f i (r + c ij t+l)-f l (r,t) = ^'°\f(r,t)), (2.20) 
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with the nonlinear collision operator fl\ ' (/) defined in ( |A.(j| ). This is the 
standard form of the nonlinear lattice Boltzmann approximation (mean field 
theory) for lattice gas automata fl[]. The way in which this equation has been 
derived closely parallels the derivation of the nonlinear Boltzmann equation 
for hard spheres from the corresponding BBGKY hierarchy ||, where one 
replaces the pair distribution function in the first hierarchy equation by a 
product of single particle distribution functions, or equivalently, one neglects 
the two point correlation function. 



2.4 Generalized Boltzmann and Rings 

In first approximation the terms linear in Gki are included in the kinetic 
equations, but higher order terms are neglected. The nonlinear Boltzmann 
equation fl2.20|) is extended with a term linear in G, yielding the generalized 
Boltzmann equation, 

f t (r + Q, t) - f^r, t) = ^'°\f(r, t)) + £ ^ k f(f(r, t))G H (r, r, t). (2.21) 

It describes corrections to the mean field equation ( |2.20|) , caused by the 
correlations Gij between the colliding particles. These correlations can be 
calculated from the lowest approximation to the second hierarchy equation. 
The derivation is given in appendix A. The result is the so-called ring kinetic 
equation, 

G y (r + a, r' + Cj ,t + l)= u ijM {f, f)G ki {r, r', t) + 5{r, r')B l3 (r, t), (2.22) 

where u is the two particle collision operator, given by 

wwnCf./O = {l + flMa^tM + fl 11,1 '^*))}^ (2.23) 

with ilf'' 1 \f{r, t)) defined in (|A.6|) . The on-node source term Bij(r,t) is 
defined by (|A.10p and ( |A.14[ ) and depends on both fi(r,t) and Gij{r,r,t). 

The ring equation ( |2.22| ) in combination with the generalized Boltzmann 
equation (|2.21| ) form a closed set of equations for the functions {fi{r,t), 
Gij(,r,r' ,t)}. These approximate kinetic equations ( [2.21 ) and ( j2.22j ) still 
obey the global (standard and spurious) conservation laws, since the orthog- 
onality condition with respect to the collisional invariants <2j, 

^ a^B l3 = 0, (2.24) 

ij 



8 



is satisfied for By, as defined by (|A.10|) and (|A.14 ). We finally note that 
^Hj' (/) i s the linearized version of the Boltzmann collision operator, ob- 
tained by Taylor expansion, in powers of 5fi, of the coefficient f^ 1 ' ^/) in 
([A. 6 ) around some arbitrary distribution fi = fi(r,t), 



^'°\f + sf) = 0?%) + £ nS ,1} (/)*/i + • • • > ( 2 - 25 ) 

j 

We discuss the assumptions and the relevance of the generalized Boltz- 
mann and ring kinetic equation for LGA's that satisfy or violate semi-detailed 
balance. The generalized Boltzmann equation and ring equation also describe 
non-equilibrium phenomena, such as transport coefficients. Application of 
the Chapman-Enskog method to the system of equations ( p.21|) and ( |2.22| ) 
will yield its normal solution and the transport coefficients. For the special 
case of detailed balance models the corrections to the Boltzmann transport 
coefficients resulting from our equations have already been calculated JTOj], 



and the agreement with simulation results is quite good ||11|| . For LGA's 
violating detailed balance the numerical calculation of transport coefficients 
beyond the Boltzmann approximation has become feasible with the present 
theory, but a detailed analysis has yet to be performed. 

Moreover, the present ring kinetic theory describes the algebraic long 
time tails of the velocity correlation and other current correlation functions, 
whereas the Boltzmann equation predicts only exponential decay. Also here 
the analytic predictions |12|] are in excellent agreement with the result of 
simulated long time tails for LGA's satisfying detailed balance [|13| . No long 
time tail measurements nor theoretical calculations have been reported in 
the literature for LGA's violating semi-detailed balance. 

The cluster expansions of the first two hierarchy equations also generate 
terms of degree b in the pair functions. In principle these terms might be 



included in ( 2.21 ) and ( |2.22|) , but a numerical solution does not seem feasible 



at present. Moreover, such extensions seem unnecessary, to judge from the 



good agreement between theory and simulation results reported above |TT 
and in the following sections. 

The equation of motion ( |2.22| ) for the equal time correlation function 
Gij(r,r' ,t) can also be derived from the Boltzmann-Langevin equation, ob- 
tained by adding a stochastic source (noise) to the nonlinear Boltzmann 
equation. In the weak noise limit the equation (|2.22| ) for Gij(r, r', t) follows, 
where the source term is the covariance of the stochastic noise source 
[|14] , [5|. The form of this source is not specified in such a phenomenological 
approach without further input, whereas the derivation here provides the 
explicit representation of Bij(r,t) as a function of fi(r,t) and Gij(r,r,t). 
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The equations above were derived under the assumption that triplet- and 
higher correlations, as well as products of pair correlations are small. The 
justification of these assumptions can only be given a posteriori by comparing 
the theoretical predictions with results from computer simulations. This will 
be done in the remainder of this paper for the special case of equilibrium 
correlations. 



3 Equilibrium correlations 

The purpose of the present section is to study the approximate equations 
for the equilibrium correlation functions, and obtain solutions under the as- 
sumption that a spatially uniform equilibrium state exists. The single particle 
distribution function will be denoted by fi(r, oo) = /j and the pair correla- 
tions by Gij(r, r', oo) = Gij(r — r'). As the on- node correlations play a 
special role, it is convenient to denote them simply by Gij = Gij(r, r, oo). 

For LGA's satisfying semi-detailed balance it is well known that the equi- 
librium distributions are completely factorized [jl], [| . Therefore the solutions 
should reduce to Gij(r — r') = S(r, r'^b^g^ This is indeed the will be 

briefly explained at the end of this section. In LGA's violating semi-detailed 
balance, Eq. ( |2.22| ) gives a quantitatively correct prediction for the on- and 
off-node correlations, extensively measured in the literature |3|, [|, [| [J . This 
is the main subject of this and the next sections. 

The kinetic equation (|2.21| ) for f± simplifies to 



^' U) +E^^ = 0. (3.1) 

k<e 

The ring equation ( |2.22| ) for the pair correlations can also be simplified. After 
introduction of the Fourier transform, 

Ee- <9T G«(r), (3.2) 
r 

it takes the form 

Gij(q) - s ij (g)a>y )JW G w (g) = s ij (q)B ij (3.3) 
or in matrix notation 

{l-s(q)u J }G(q) = s(q)B, (3.4) 

where G(q), G, and B are 6 2 -dimensional vectors with components Gij(q), 
etc., and uj, s(q), and 1 are b 2 x b 2 matrices with elements u)%j,kti s ij,ke(9) = 
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Sij(q)5ikSje, and lij t u = Sik$ji- The pair streaming operator is denned as 

Sij(q) = exp[-iq • (Ci - Cj)]. (3.5) 
In the spatially uniform equilibrium state the quantity B^ simplifies to, 

Bij = + ^2 ^ijlki^ki + XX ^ ~ ^ij^Gu- (3-6) 

k<l kl 

Note that due to (|3.1|) the contribution B^ in (|A.10| ) vanishes identically. 

In appendix B we discuss the technicalities of obtaining a well-defined 
closed equation ( |B.12| ) for the on-node pair correlation function G^ or more 
conveniently for the excess correlation function Cy = Gij — S^gi. This equa- 
tion reads 

C = R{Q^ + (1 - u)}C + M2 (2 ' 0) + J, (3.7) 

where the ring operator R and the term J are given by ( B.10|) and (B.13| ) 
respectively. Most technical difficulties are related with finite size effects that 
give deviations for small system sizes V = L d . In the thermodynamic limit 
(V — > oo) many complications disappear: the term J vanishes, and the ring 
operator becomes 



R = lim 



1 1 f ^ f dq 



V^oo 



qeiBZ 



V nTT^r, 1 — s(q)u jibz (2-n) d 1 — s(q)u 



where summation and integration are restricted to the first Brillouin zone 
(1BZ) of the reciprocal lattice C* and where vq denotes the volume of a 
unit cell on the lattice £ (e.g. v — 1 for the square and Vq = \\f^> for the 
triangular lattice). 

Once the on-node correlations are known from the solution (|3.7| ), then 



Gij(q) can be calculated from (|3.4j ), and the full r-dependent correlation 
function Gij(r) follows by inverse Fourier transformation of Gij(q), 

^(r)4E^«(9)- ( 3 - 9 ) 
v Q 

In equilibrium the pre- and postcollision correlations are related by 

G? j {r) = G ij (r + c i -c j ), (3.10) 

because of the relation n,(r + Cj, t + 1) — n*(r, t) between pre- and postcol- 
lision occupation numbers. 

Before concluding this section we make some further comments on the 
ring equation ( |3.7|) . The term RQ( 2 ' ^ on the right hand side yields the simple 
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ring approximation to the (precollision) correlation function C. The quantity 
^ij'°\f°), defined in ( |A.15| ), constitutes the lowest order approximation to 
the postcollision correlations. It represents the on-node postcollision corre- 
lations, created by a single collision from a completely factorized precollision 
state. This quantity was calculated earlier in Ref. and agreement within 
10% was observed in computer simulations on a triangular lattice gas violat- 
ing semi-detailed balance. The present theory enables one to extend these 
calculations to the precollision correlations, which are built up by collective 
phenomena, acting on large spatial and temporal scales. The ring opera- 
tor ( p.lO|) propagates these two particle correlations to other nodes through 
uncorrelated single particle motion (Boltzmann propagators), and recollects 
this information on a single node. 

The term R{Q( 2 ' 2 ^ + (1 — uj)}C on the right hand side of (3/7) represents 
the repeated ring collisions. Numerical evaluation of ( |3.7| ) in the next sec- 
tion will show that the repeated rings represent only small corrections to the 
equilibrium correlations, as calculated from the dominant simple ring approx- 
imation, RQ( 2 ' 2 \ Nonvanishing equilibrium velocity and position correlations 
are only found in LGA's that violate the condition ( |2.4| ) of semi-detailed bal- 
ance. We will return to such LGA's in the next section. 

In the remainder of this section we show that the excess correlation func- 
tions Cij are vanishing in LGA's that satisfy the semi-detailed balance condi- 
tion. If mass and possibly momentum is conserved in a LGA, the equilibrium 
solution is fi = p/b with p the average number of particles per node. This can 
be verified from the nonlinear Boltzmann equation, showing (/°) = 0. 

The factorized node distribution (|2.19|) obeys the relation, F (s) = F (a), 
where s and a are respectively pre- and postcollision occupation numbers. 
This is so because F (s) depends only on the conserved number of particles, 
J2i s i( r ) — J2i °"i( r ) on a single node. The relation F (s) = F (a) also implies 
that fig' 0) (/°) = in ( |Og ) for semi-detailed balance models. This can be 
verified by using in the first and second term on the right hand side of ( |A.15| ) 
respectively the semi-detailed balance condition ( |2.4| ) and the normalization 
condition ( |2.2j ). Consequently, the ring equation ( B.12|) gives Cy = for 
semi-detailed balance models. There are no on- or off-node velocity correla- 
tions in the equilibrium state of LGA's that obey the semi-detailed balance 
condition, i.e. Gij(r,r f ) = 5ij5(r,r')f®(l — /P), in full agreement with the 
Gibbs distribution. 
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4 Numerical evaluation 



So far, we have derived an equation to determine static pair correlations. 
This section deals with the numerical evaluation of the theory from a more 
operational point of view. We discuss an iterative scheme for finding the 
equilibrium solution numerically, and make some remarks about its unique- 
ness. 

The stationary solution {/, G} must satisfy the generalized Boltzmann 
equation and the ring equation (|3.7|) . In single speed LGA's violat- 



ing semi-detailed balance lattice symmetries require that the single particle 
distribution function is given by f® = p/b. To calculate the on- node correla- 
tion function from (577) we first evaluate f2( 2 ' ) and f^ 2,2 ) for ff = p/b, and 



then solve the linear matrix equation ( [B.9| ) to obtain the excess correlation 
function Cy. 

However, in multi-speed models (for instance, with additional rest parti- 
cles) that violate semi-detailed balance the stationary single particle distri- 
bution generally differs from p/b. How can one numerically determine the 
equilibrium solution in the general case? We have found it efficient to use 
the following scheme, which amounts to self-consistent iteration of (|3.1|) and 



1. Specify an accuracy e C 1. 

2. Find the stationary distribution ff in Boltzmann approximation, i.e. 
with all pair correlations Gij set equal to zero. Starting from = 
p/b we repeat the iterative step, f(n + 1) = f(n) + f2( 1,0 )(/(n)), until 

< e. Note that f^ 1 ' -* must be re-calculated at every iteration 

step. 

3. Calculate the (unique) solution of the ring equation (|3.7| ) for given 
fixed /P. 



4. Find the stationary solution ff of (|3.1|) for given fixed Gy as found in 



step (3). This can be done by repeating the iterative step f(n + 1) = 
f{n) + I(/(n)), with / = tt( 2 > *> + tt^G, until |/| < e. Again I must 
be re-calculated at every iteration step. 

5. Repeat steps (3) and (4) until the scheme has converged. 

Step (3) in this scheme requires some explanation. The ring operator R is 
calculated using the spectral decomposition ( |B.2|) . This requires a determi- 
nation of the q-dependent eigenvalues X a (q) and the complete bi-orthogonal 
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set of left- and right eigenfunctions {Xa(q), Xa(q)} of the two particle prop- 
agator 7(q) = s(q)u. Once R is known, we proceed to calculate Gij or 
for given /? by splitting the ring equation ( |3.7|) into an inhomogeneous part 
that depends only on ff, and a part that depends only on the excess pair 



correlation functions. We can write (|3.7|) as 



C = MC + K. (4.1) 

with M = R(ttW + 1 - oj) and K = Rtt (2 >V + J. Since C tj = is sym- 
metric, the only independent elements of are those with i < j. These can 
conveniently be taken together as a vector in a \b{b — l)-dimensional space, 
spanned by the possible pairs (ij) with i < j. It is therefore appropriate to 
interpret (3~T) as a vector equation, where C and K are \b(b — l)-dimensional 



vectors and M is a matrix of equal dimensionality. Once K and M are known, 
C = (1 — M) K can readily be calculated, provided that det(l — M) 7^ 0. 
As noted earlier, the term J vanishes in the thermodynamic limit. 

At this point we have completed the detailed description of our scheme 
for finding the stationary solution to the coupled time evolution equations 
for fi(t) and Gij(t). Once the equilibrium values of ff and Gy have been 
found, Gij(q) and G(r) can be calculated from ( |3.4|) and (|3^ ). One can now 
ask whether the solution that is obtained with this scheme is unique. The 
stationary Boltzmann distribution ff, obtained in step (2), is by definition 
in the basin of attraction of the fixed point {/°, G} found using our itera- 
tive scheme. In principle it can not be excluded that there are other stable 
fixed points that correspond to physically acceptable equilibrium states as 
well, since the f2-matrices are nonlinear functions of / = However, 
we can systematically search for other fixed points, since for a given / the 
corresponding Gij(f) is uniquely determined by the ring equation. Then by 
investigating Ii(f,G(f)) = S^ 1 ' ) + Q^' 2 >G, as a function of / we can locate 
the fixed points, given by the condition Ii{f, G(f)) = 0. 

It is important to stress the difference between two alternative ways of 
finding a spatially homogeneous stationary solution to the time evolution 
equations, Q2.2ip and fl2.22| ). The first and most efficient method uses the 



iterative scheme discussed in this section, which can be interpreted as a way 
to find the fixed point to a mapping {/("), G^} -> {f {n+1 \ G^}. A second 
method is to specialize Q2.2H) and (|2.22|) to the spatially homogeneous case, 
which yields another mapping, {/(£), G(r, t)} — > {f(t + l),G(r,t + 1)}. If 
there is more than one stationary solution, then the two schemes will not 
necessarily converge to the same fixed point {/, G}. 
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5 Applications 

5.1 Interacting random walkers 

The ring kinetic theory for approximate calculation of pair correlations pre- 
sented in the previous sections assumes that all higher order correlations are 
negligible. It is important to make a quantitative comparison with computer 
simulation results to establish the accuracy of the results obtained from a 
numerical evaluation of the theory. 

We first consider a model of interacting random walkers on a line. An 
isolated random walker executes a persistent random walk, while jumping to 
nearest neighbor sites or resting on the same site, i.e. it has three allowed 
velocity states, Cj = {+1,0,-1} with transition probabilities depending on 
the previous jump, as defined in the left diagram of Fig. |l|. 

The interaction between the walkers is strictly local, and essentially de- 
termined by the Fermi exclusion rule for different velocity channels, i.e. 
rii(r,t) = {0,1}. The maximum number of walkers on a site is therefore 
three. The transition probabilities from a two-particle in-state s(r) to an 
out-state a(r) are defined in the right diagram of Fig.l. There exists only 
a single three-particle state, which remains unchanged under interaction on 
account of the Fermi exclusion rule. 

The above model can be described conveniently in terms of a one-dimen- 
sional diffusive LGA, where at every node one- and two-particle transitions 
occur with probabilities defined in Fig.l. In principle there are twelve inde- 
pendent transition probabilities A sa ^ 0. This number is reduced to six if we 
require that the collision rules be invariant under reflection (q — > — Cj). One 
may further impose self- duality, i.e. invariance under exchange of particles 
and holes (n^ — > n = 1 — n^) . This would reduce the number of independent 
transition probabilities to three, so that a — a', (3 — /?', 7 = 7'. In the present 
model the total number of particles is conserved, but total momentum is not, 
and there are no staggered or other spurious invariants. Consequently the 
model has only a single slow (diffusive) mode. 

To illustrate some of the analytical results of the previous sections, and 
to discuss the conditions of semi-detailed balance, it is of interest to write 
out the microdynamic equation ( |2.6| ) for the occupation numbers n + (r,t), 
no{r,t), and n_(r, t), referring respectively to the velocity channels {q = 
+ 1,0,-1}. Let n*(r,t) = riiir + q,£ + 1) be the postcollision occupation 
number, then 

n* + = n + + (d n n + — (3 n + n )n_ + 7 {nJn + — n + n_)n + 
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= n o + ft {n+n- — n_n + )no — ct r n n + n^ — ai rion + n- 

+a' r n + n_n + ajn + n_n — ft' (n + n_ + n_n + )n , (5.1) 

and a similar relation for n*_. The set A sa = {a r , a.i, a' r , a' h ft, ft', 7, 7'} rep- 
resents the Boolean realizations, A sa (r,t), discussed below (|2.6| ), with ex- 
pectation values A sa = {a r , an, a' r , ot[, ft, ft' , 7, 7'}. The subscripts {r, 1} refer 
to transitions to states with a particle moving to the right and to the left, 
respectively. 

Consider first the mean field or Boltzmann approximation. According to 
( |2.20| ) the nonlinear Boltzmann equation for fi(r,t) is obtained from (|5 . 1| ) 
by replacing n,(r, t) by fi(r,t) and the Boolean variables a, etc by their 
expectation values a, etc.. In equilibrium the distribution function ff is 
characterized by the average occupation numbers /o and / + = /_ with av- 
erage density p = fo + 2/ + . At a given density p the nonlinear Boltzmann 
equation for the stationary state reduces to a quadratic equation for / . 

Inspection of Q5.1D shows that the semi-detailed balance condition Q2.5Q 
is only satisfied if 

««' = pp', (5.2) 

in which case the stationary distribution is given by z /z + = (3 /a = a'/ ft', 
with Zi = /i/(l — fi)- We recall that the stationary distribution is a Gibbs 
distribution if the phase space density depends only on the conserved quanti- 
ties, so that fi=fo = f+ = f~ = lp is independent of the velocity channel. 
Therefore detailed balance with respect to the Gibbs distribution imposes 
the condition, 

a = ft, a' = ft'. (5.3) 

With the small set of allowed states s(r) the condition for semi-detailed 
balance reduces to that for detailed balance [|]. 

If the above detailed balance conditions are violated, there exist correla- 
tions between the occupation numbers of the different channels (r, Cj), even 
in the equilibrium state. The critical quantity, that determines whether equi- 
librium correlations are non- vanishing, is matrix element Q^'°\f°), defined 
in ( |A.15| ). It gives for the present model, 

4 2 L°V) = 2(ft'f J + -a'f + f + )U 

n%°\f) = n { - \f ) = (a'f + J -ft'f J + )f + . (5.4) 

If the transition probabilities satisfy the detailed balance conditions, then 
^if°\f°) is vanishing for Vi, j. All excess on- and off-node correlations van- 
ish, and the stationary distribution is simply /o = /+ = /_ = |p. In case 
the detailed balance conditions are violated, the elements of Q^ 2 '°'(f°) differ 
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from zero, and there exist non-vanishing on- and off-node correlations in the 
equilibrium state. 

The method developed in the present paper enables one to calculate the 
equilibrium distributions {fi,Gij} numerically, from which the off-node cor- 



relations Gij(r) can be constructed with the help of ( |3.9|) . The postcollision 
correlations Cr*-(r) are then given by ( |3.10| ). Note that for the 3-bit model 
G*_ + = Cry (2) and Gq + = Cr +(1), where r = 1 and r = 2 denote the nearest 
and next nearest neighbor sites. The relative importance of the pre- and 
postcollision pair correlations Cry and GL is best measured when they are 
normalized by the single channel fluctuations, gi = ((<5nj) 2 ), yielding the 
covariances^ 

Cav(M) = -= Cov*(i,j) = -^. (5.5) 

Using the iterative scheme of section 4 we have calculated these quantities 
for different choices of reduced density f = p/3, lattice size L, and transition 
probabilities a, j3, 7 (for self-dual models only), and compared the theoretical 
results with computer simulations. 

In Fig. ^|a values of the pre- and postcollision on-node correlations, ob- 
tained both from the numerical evaluation of the theory (lines) and from 
computer simulations (symbols with error bars), are plotted as a function of 
the system size L. For small values of L there are strong finite size effects, 
which are quantitatively very well predicted by the theory. 

It should be stressed that in a closed (microcanonical) system, with finite 
L and fixed N, correlations are always present between occupation numbers, 
even in models satisfying semi-detailed balance. This can be seen as follows: 
for a semi-detailed balance model with V nodes and N = bfV particles 
we have if (r*,Cj) 7^ (r',Cj) the following relation, (ni(r)rij(r')) = f[(N — 
l)/(bV - 1)] < f 2 , so that Gij(r) = -0(1/V) < 0, independent of r. 

Fig. 0b shows the typical dependence on the reduced density / = p/3; 
the correlation functions are symmetric around / = 0.5, as a consequence of 
the imposed self-duality. In Fig. ||c the on-node correlations are plotted as a 
function of a, for fixed (3 = 0.33 and 7 = 0.5 at reduced density / = 0.5 and 
L = 128. Note that when a — [3 all correlations vanish: the model satisfies 
(semi)-detailed balance. 

Since there is no fundamental difference between on- and off-node correla- 
tions, and since in addition the L-dependence of the correlations is well pre- 
dicted by the theory, we expect that our theory is capable of giving good pre- 

2 In Ref. [|J the correlations Gij were normalized by fafj rather than by y/gjg], as is 
done in this paper. Only with the latter choice the self-dual symmetry is preserved. 
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dictions for off- node correlations Gij(r) as well. This expectation is justified 
by Fig. [3]a,b, where for fixed / = 0.5, L = 128 and two different sets of transi- 
tion probabilities, the density- density correlation function G(r) = J2ijGij(r) 
is plotted versus r. Note that G(r) is negative at large r, due to the finite 
size effects. 

To study how the system approaches equilibrium in a spatially uniform 
non-equilibrium state we have solved the time-dependent evolution equations 
( |2.21|) and (|2.22|) . Figure |]a shows how the typical relaxation time increases 
with L. Note that for t < L/2 the evolution is independent of L, which 
can be explained by the observation that the shortest time in which two 
opposite moving particles can meet through the periodic boundary conditions 
is t = L/2; at earlier times there is no distinction between a finite and an 
infinite system. When L = 1024 the system still has not reached equilibrium 
after 1000 time steps. 

For long times the approach to equilibrium is algebraic, Gy(t) — Gy(oo) = 
AGij(t) ~ t~ a , with an exponent a that depends on the dimensionality and 
the type of collisional invariants of the model. Fig. shows that a = | in 
the case of the interacting random walkers, where the number of particles is 
the only collisional invariant. In a separate publication we will analyze the 
exponents and amplitudes for these algebraic tails in detail. 

5.2 Two-dimensional fluid-type model 

In a previous publication Q on-node correlations have been studied in a 7- 
bit two-dimensional fluid-type LGA on a triangular lattice, which violates 
(semi)-detailed balance. Each node can contain a rest particle, and up to six 
moving particles with velocities corresponding to nearest neighbor vectors. 
Even after imposing that the collision rules conserve the lattice symmetries, 
the model still has 20 independent transition probabilities. 

Special choices of these probabilities, referred to as set #n, have been 
studied in Ref. |Q, where on-node correlations were measured from computer 
simulations. For certain choices of transition probabilities (e.g. set #30) 
the correlations did not reach an equilibrium value, but kept growing in 
time. This behavior could later be assigned to the existence of unstable 
sound modes, driving a phase separation However, there were other 

choices of transition probabilities (e.g. set #13) that did violate detailed 
balance while all modes remained stable. In these stationary state 

was reached after a few hundred time steps. In that paper only the on- 
node postcollision correlations could be calculated in lowest approximation, 
i.e. C*j = Qij'°\f°) (see discussion below ( |3.10| )), but a quantitative theory 



18 



for the precollision correlations in equilibrium was entirely lacking Q. The 
present paper provides the missing theory. 

To test our theory for this non- detailed balance LGA we have chosen the 
set #13 of transition probabilities, and evaluated the ring kinetic theory The 
theoretical predictions for the correlation functions are compared with new 
computer simulations of higher statistical accuracy, taking several hours of 
CPU time on a RISC workstation for each data point. In basic equilibrium, 
where the total momentum P vanishes, the only independent elements of 
Cij are (ij) = {(01), (12), (13), (14)}; all other elements are related to these 
four by lattice symmetries. As set #13 corresponds to a self-dual model, the 
correlations are again symmetric around / = p/7 = 0.5. Fig. |5]a-d shows 
the simulation results (symbols) compared with the ring kinetic theory for 
precollision (solid line) and postcollision (dashed line) correlations. From 
this figure it is clear that the self-consistent ring kinetic theory agrees very 
well with the simulations. For comparison we have also plotted the functions 
^ij'°\f°) (dotted line) which represent the postcollision correlations created 
by a single collision step acting on a factorized precollision state with average 
occupations f° obtained from the nonlinear Boltzmann equation (see discus- 
sion below eq. |2.25|) . In an earlier paper j|] we obtained C*j = f^ 2 '°)(/°) as 
a simple estimate for the postcollision correlations. 



6 Conclusion 

We summarize the conclusions, consequences and possible applications as a 
set of comments. 

1) In the present paper, for the first time a theory is presented to calculate 
distribution and correlation functions in lattice gas automata (LGA's) that 
violate semi-detailed balance. The theory of standard LGA's satisfying semi- 
detailed balance can be recovered as a special case. 

In section 2, starting from an open BBGKY hierarchy of coupled time 
evolution equations for the n-particle distribution functions, we obtain clo- 
sure by neglecting correlations of third and higher order. We are left with two 
coupled time evolution equations for the single particle distribution function, 
fi(r,t), and the pair correlation function, Gij(r,r' ,t). In section 3, to gain 
insight in the equilibrium state for LGA's violating semi-detailed balance, 
we specialize our theory to a spatially uniform stationary state, fi(r,t) = 
and Gij(r, r', t) = Gij(r — r'). For given fa, the on-node correlations are 
uniquely determined by a linear equation: the so-called ring equation (|3.7|) . 

2) The ring operator R defined in (|3.8|) is identical to the one introduced 



for calculating long time tails of time correlation functions [12, 16] and cor 
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rections to Boltzmann transport coefficients |10| in LGA's satisfying semi- 
detailed balance. One can show that linearization of the theory in this paper 
around the Gibbsian equilibrium state for standard LGA's reproduces the re- 
sults of Refs. [JT2], [0], |IIJ . Using the spectral decomposition method discussed 
in appendix B - rather than the method discussed in Ref. 0] - in combi- 



nation with a Gaussian integration scheme for performing the q-integral in 



(3.8), the ring operator can efficiently be evaluated numerically for large 
system sizes. 

3) In section 4 we provide a numerical scheme for finding, by self-consistent 
iteration, a stationary solution {fi,Gij} that is 'close' to the solution of the 
nonlinear Boltzmann equation (where Gij = 0). Is the equilibrium solution 
found by this scheme unique? We note that even though the ring equation 
uniquely determines Gy for given f i: it is highly nonlinear in Therefore 
it can in principle not be excluded that more than one spatially uniform 
stationary solution, or equilibrium state, exists. 

4) It is well known that the violation of semi-detailed balance in LGA's may 
lead to spatial instability of the uniform state and subsequent pattern forma- 



tion UTTf . Computer simulation studies of a two-dimensional LGA violating 



semi-detailed balance, where the presence of unstable sound modes leads to 
phase separation, have indicated that a finite system with periodic bound- 
ary conditions randomly 'chooses' between different spatially inhomogeneous 
stationary states, corresponding to limit cycles with a 'continuously' broken 
symmetry |15[ . 



5) To test our theory we have applied it to two different LGA's, both violat- 
ing (semi)-detailed balance: a one- dimensional model of interacting random 
walkers, and a two-dimensional fluid-type LGA defined on a triangular lat- 
tice. For the first model we have compared theory with simulations for var- 
ious densities, transition probabilities and system sizes. In all cases there is 
remarkable agreement. Even the finite size effects occurring for small system 
sizes are well predicted by the theory. The small deviations that occur in a 
few cases are probably due to the fact that third and higher order correlations 
are systematically neglected. 

6) The second model was already considered in an earlier publication [|]]. 
In that paper a simple theoretical estimate was given for the postcollision 
on- node correlations, C*j = fiy , in reasonable agreement with computer 
simulations. However, theoretical predictions for the precollision on-node 
correlations, and for off-node correlations were totally lacking. The present 
paper provides the missing theory, and also shows the excellent agreement 
between theory and simulations for this two-dimensional fluid model. 

7) As t — > oo a stationary state is approached. Due to the existence of local 



20 



conservation laws the approach to this stationary state is algebraic, ~ t~ a , 
with an exponent a that depends on the dimensionality and on the type 
of conservation laws. There is an intimate connection between the algebraic 
tails in the approach towards the correlated equilibrium state mentioned here, 
and the well-known long time tails in the decay of the velocity autocorrelation 
function [JT^J, [Tj|. More generally, several ideas from the field of classical 
kinetic theory can be applied to the equations derived in this paper, which 
will be discussed in a separate publication. 

8) In Ref. |14| the Boltzmann-Langevin equation for LGA's was discussed. 
This phenomenological equation requires specification of the correlation func- 
tion of the random noise sources. To what extent does the fluctuation- 
dissipation theorem enable us to relate the noise strength to the phenomeno- 
logical damping or transport coefficients? The fluctuation-dissipation the- 
orem requires the stationary distribution and correlation functions as in- 
put. If one does not know this distribution function a priori the fluctuation- 
dissipation theorem is an empty statement. For LGA's that violate (semi)- 
detailed balance, one may therefore also summarize the implications of the 
previous discussion by stating that the standard form of the fluctuation- 
dissipation theorem for Langevin's equation, relating the noise strengths to 
the dissipation coefficients, is no longer valid, because the stationary distri- 
bution differs from the Gibbs distribution. 

9) It is appropriate to give some comments on the validity of the generalized 
Boltzmann equation and ring equation derived in this paper. In real fluids 
the validity of the cluster expansion is restricted to dilute systems, where the 
density acts as a small expansion parameter. Here there is no such small pa- 
rameter and the approximations are less controlled. To understand the great 
success of the present mean-field theory, it is useful to compare the lattice 
Boltzmann equation with the revised Enskog theory for hard sphere fluids 
[ |T9|j . Both theories describe non-equilibrium phenomena reasonably well over 
the whole density range. The reason for this seems to be that the collision 
term for the Revised Enskog Theory (resp. Lattice Boltzmann equation) con- 
tains the static correlations imposed by the hard core exclusion (resp. Fermi 
exclusion), which are the origin for the strong density dependence of the the- 
ory. The present ring kinetic theory equation also contains - through the 
expansion coefficients Q - all static correlations between colliding particles 
imposed by the Fermi exclusion. This might explain the success in predicting 
the density dependence of the correlation functions. 
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Appendix A 



In this appendix the cluster expansion of the hierarchy equations is derived. 
We start by substituting n» = /» + 5n,i into identity Q2.13Q and obtain, 



(/< + 5n^(l - ./• - Sntf-i = /; (1 - f 3 ) 



\ l — Sj 



9j 



(A.l) 



with gj = fj(l — fj). With the help of this identity we can expand the single 
node distribution p(s, r,t), defined in 



p{s, r, t) = F(s, r, t) ll + £ S J^^S. Gh , {r .r.t) + . . .} . (A. 2) 



9k(r,t)gi(r,t) 



where r, t) is the completely factorized single node distribution function 
defined in ( |2.19| ) . Similarly we can expand the two-node distribution function 
in ( frig ) as 



p ( J (s, r, s', r', t) - p(s, r, t)p(s', r , t) = 

F(s, r, t)F(s', r>, t) £ ^j^ G M (r, r',t) + .... (A.3) 



# fe (r,t) gi(r',t) 



Using the normalization ( |2.2| ) of A scr we rewrite the time evolution equation 
( |2.11|) for the single particle distribution function fi(r,t) as, 



fi(r + Ci, t + 1) - /i(r, t) = Yj( a i ~ Si)A sa p(s, r, t). 
Insertion of ( A.2|) yields 



(A.4) 



fi{r + Ci ,t + 1)- f t {r,t) ~ ^ (i ' Uj (/) +EnS J (/)G M (r,r ) t). (A.5) 
The expansion coefficients f2 are given by 

«S 1,0) (/) = E(^-^K^) 



5sj 



so- 9j 



(A.6) 



where k ^ £ and -F(s) = F(s, r, t) is the completely factorized single node 
distribution function defined in (j2.19|) . Note that the f2's depend on r and t 
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through f(r, t). They are generated by the Taylor expansion of the nonlinear 
Boltzmann collision operator, fq'°\f + 5f), in powers of 5f. 

To derive an equation of motion for Gij(r,r',t), defined in ( 2.16Q , we 
start from 

djir + Ci^r' + Cj ,t + 1) = f§\r + c i ,r' + c j ,t+ 1) - (A.7) 
f i (r + c i ,t + l)f j (r , + c j ,t+l) 

and use ( 2.11 ) together with ( |2.12 ) to obtain 

G? y (r + (k, r' + Cj , t + 1) = K i3 {r, r', t) + 5{r, r')B i3 (r, t), (A.8) 

with 

Kij(r, r', t) = J2 E SviSv'jAsaAs'a, {p (2) (s, r, s', r', t) - p(s, r, t)p(s', r', t)} . 

(A.9) 

and the on-node source term, 

B i3 {r,t) =BW(r 1 t)+B§\r,t)+B§\r 1 t) 1 (A.10) 
consisting of three parts, 
B§\r,t) = ^(5difo,4p(s,r,i) 

= Gij(r, r, t) + ^(SaMj - 8sidsj)A scr p(s) 

StJ 

B§ (r,t) = SvM'jAsaA^, [pV ( s , r , s', r, t) - p(s, r, t)p(s', r, t)} 

SIT s'o-l 

B<?\r,t) = -J2Y, 6 ^ 6 °j A s« A s'Ms,r,t)p(s',r,t). (A.ll) 

su s 'a' 

The cluster expansion of K i3 follows from ( |A.3| ) with the result 

K i3 (r, r', t) ~ to ijM (f, f')G ke (r, r', t), (A. 12) 

where the pair collision operator is defined as 

^ jM (f,f) = E SaMjAsaA^FisMs') 5 -^ 

sas'a' ' 9k9e 

= {l + n^(/(r, + n^(/(r / ,t))}^. (A.13) 
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To identify f^ 1 ' 1 ^/) with the symbol introduced in ( |A.6| ) one needs to verify 
that the term containing SsiSsj equals Sij . For the three parts of Bij(r, t) we 
obtain 

B§\r,t) ~ G 1J (r 1 r,i) + nf(/) + ^ng(/)GH(r 1 r 1 t) 

k<e 

Bg\r,t) ~ -w y , w (/,/)G w (r,r,t) (A.14) 
if (r,t) - -^ 1 ' 0) (/)^° ) (/) 

- E fa™ (/)«&?(/) + «T (/)«&?(/)} G«(r,r,t). 

We have introduced two more expansion coefficients, where k 7^ £, 

= E^-M^W^. (A.15) 



Appendix B 

In this appendix we discuss the formal solution of the ring equation 

G(q) = 1 s(q)B. (B.l) 
1 — s{q)uj 

Using a bi-orthogonal set of right and left eigenvectors, Xa(o) an d Xa(o), 
and eigenvalues X a (q) of the two-particle propagator 7(g) = s(q)u, we can 
write (1 — s(q)u)^ 1 as a spectral decomposition, 

T 3^ = i:ii»M) T 3^y(^)l, (B.2) 

provided that X a {q) 7^ 1. There is however a complication if, for one or more 
reciprocal lattice vectors q n belonging to a set {<?„}, there are eigenmodes 
a for which \ a (q n ) = 1. The operator (1 — s(q n o;)) then has a null space, 
spanned by the eigenvectors for which X a (q n ) = 1. The set of points {<?„} 
contains the origin, corresponding to the standard conservation laws, and the 
centers of the facets of the Wigner-Seitz cell, corresponding to the staggered 



invariants [20] . 



We can now formulate a necessary solubility condition, which requires 
that the inhomogeneous term in (|3.4|) be orthogonal to the null space, i.e. 



(X*(qJ \s{q n )B) = (B.3) 
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for the modes (a,n) that have X a (q n ) = 1. It can be verified that this 
condition is indeed satisfied as a consequence of ( [2.24 ), i.e. the local conser- 



vation laws, since both for the standard conserved quantities TC and for the 
staggered invariants we have s(q n ) = 1. 

After establishing that ( |3.4f ) is soluble, we observe that its solution is not 
unique in the set of points {q n }, where it is given by 

G(q n ) = Qn^—B + ^SnaXaiQn), (B.4) 

where the relation s(q n ) = 1 has been used. The coefficients £ na are undeter- 
mined, and the prime on the sum indicates that the summation is restricted 
to eigenmodes a for which \ a = 1. The projector Q n , defined as 

Q n = 1 - P n = 1 - Y,'\Xa(q n ))(Xa(q n )\, (B.5) 

a 

is introduced to make explicit that the inverse of (1 — s(q n )uj) is only defined 
in the orthogonal complement of the null space of (l — s(q n )u), where X a ^ 1. 
The second term lies in the null space, i.e. 

PnG(q n ) = Y! S naXa(q n )- ( B -6) 

a 

The coefficients S na are in fact equal to the covariance of the fluctuations of 
the (standard and spurious) global invariants. For instance, for the standard 
invariants TC = {N,P,...}, associated with q = 0, the eigenmode x«(0) = 
a\a^ is a product of collisional invariants. Then 

r ij v 

Therefore the coefficients S na are constant in time, and can be calculated 
from the initial ensemble at t — 0. In the calculations of this paper we 
assume that the initial ensemble is prepared microcanonically, i.e. with all 
global (standard and spurious) invariants fixed, from which it follows that 

In the previous discussion we have explained the structure of the formal 
solution G(q) of (|3.4j). However, this is not yet a solution of the ring equation, 



or even a closed equation for G(q), because the inhomogeneous term Bij is 
given in terms of the on-node distribution function Gij. It is related to Gij(q) 
by 

G« = iE fi «(9). (B-8) 

v q 
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where the g-summation runs over the first Brioullin zone in the reciprocal 
lattice £*. 

Combination of (|B.8|) , QB.4Q , (|B.6|) and ([B.l|) gives the following closed 
equation for the on-node correlation matrix G, 



(B.9) 



with B given in ( |3.6| ) and the ring operator R is according to ( |B.2j ) and ( [B.4j ) 
defined as 



R 



1 

V 



S r 

qmj 



s{q)u 



■AQ) 



1 



V 



00 



(B.10) 



where s(q n ) = 1 has been used. The asterisk on the summation sign is 
defined through the second equality. 

As the diagonal part Gd of the on-node correlation function is known for 
given /?, it is convenient to introduce the excess correlation function C = 
G — G(ii and to absorb the contribution from Gd into a known inhomogeneous 
term. To do so, we observe that VlP^Gd = on account of ( |A.15 ), and we 
calculate 



AQ 



)(l-«)G*=^£*Gd. (B.11) 



R(l-cj)G d =-^2- „ 

V q 1 - S (<?V V q 

In the second equality we have used the relation s(q)Gd = Gd because 
Sa(q) = 1. Equation (|B.9|) yields then for the excess correlation function, 



C = R{n^ + (1 - uj)}C + M2 (2 ' 0) + J, 
with the finite size correction term J given by, 



(B.12) 



V 



= ^T. p n[G{q n )-G d ] 

' n 
V na 

The diagonal fluctuation is defined through, 

V na = (Xa(q n )\G d ) = Ex«Ji(9„)^(l " fj) 



(B.13) 



(B.14) 
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Note the difference between (|B.7|) with f = p/b and the present equation 
where /° is the stationary solution of the generalized Boltzmann equation 
( |2.21|) . In semi-detailed balance models f® = f, but in multispeed LGA's 
lacking semi-detailed balance /? differs from / already in mean field approx- 
imation ([2.20 ) . 

Inspection of ( |B.13| ) shows that the na-summation contains only as many 
terms as there are global invariants; so J in ( |B.12| ) is O^V" 1 ) and accounts 
for finite size effects. In the thermodynamic limit J vanishes and the ring 
equation ( |B.12| ) yields a closed equation for the excess correlation function 
C. The diagonal part Gd then no longer appears in the ring equation. 
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Figure 1: Transition probabilities in one-dimensional model of interact- 
ing random walkers with mutual exclusion in the same state {r, Cj} with 
Cj = {—1,0, 1} = {<—,«,—>}. The left and right diagrams show transitions 
between states with one and two particles, respectively. 
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Figure 2: Comparison of computer simulation results (symbols with error 
bars) and self-consistent ring kinetic theory (lines) for the model of Fig. |l[ 
On-node correlations, indicated as covariances, are plotted (a) versus system 
size L, for transition probabilities a = 0.4, /3 = 0.5, 7 = and density 
/ = 0.5; (b) versus density /, for L = 128, a = 0.5, f3 = 0.4, and 7 = 0; (c) 
versus transition probability a, for / = 0.5, L = 128, (3 = 0.33, and 7 = 0.5. 
Note the vanishing of correlations when a = j3 (detailed balance). 
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Figure 3: Off-node density- density correlation function G{r) = J2ijGij(r) 
plotted versus r, for / = 0.5, L = 256, and two different choices of transition 
probabilities: (a) a = 0.1, [3 = 0.5, 7 = 0.5, and (b) a = 0.5, (3 = 0.1, 
7 = 0.5. Simulations (symbols with error bars) compared with theoretical 
values (solid line connecting points r = 0, 1, 2, ..). 




Figure 4: Time-dependence of on-node correlations, calculated using ring 
kinetic equation (no simulations), with / = 0.5 and the transition probabili- 
ties of Fig. |2|a: (a) short time behavior for various L; (b) algebraic long time 
behavior, AG« ~ t~ 1/2 , for L = 4096. 
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Figure 5: Test of ring kinetic theory for the model of Ref. using the 
transition probabilities of set #13. Simulation data for L = 64 and / = 0.5 
(symbols with error bars), compared with old theory for on- node postcollision 
correlations of Ref. M (dotted line), and with ring kinetic theory for the on- 
node precollision (dashed line) and postcollision (solid line) correlations. 
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